Salmon VBGF model v3

Load libraries

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools","viridis","nls.multstart", "broom", "patchwork", "nimble", "coda", "boot", "tidybayes","bayesplot")

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  }

invisible(lapply(pkgs, library, character.only = T))

options(ggplot2.continuous.colour = "viridis")
theme_set(theme_light()) # check why global theme option not working when rendering
# Set path

home <- here::here()

Read data

sallaa <- readRDS(file = paste0(home,"/data/data-for-2-2/salmon-laa_2025-04-11.RData")) %>%
  #filter out non-aged individuals (~34000 individuals)
  filter(!(is.na(age.sea) & is.na(age.sm))) %>% 
  rowwise() %>% # rowwise needed for non vectorised function like sum
  # create life stage and total age columns
  mutate(age.type = case_when(age.sea == 0 | is.na(age.sea) ~ "smolt.only",
                               age.sea > 0 & age.sm > 0 ~ "both",
                               age.sm == 0 | is.na(age.sm) ~ "sea.only",
                               .default = NA),
    #mutate(age.type = if_else(age.sea == 0 | is.na(age.sea), "smolt", "sea"),
         age.tot = sum(age.sea,age.sm, na.rm = TRUE)) %>%
  ungroup()

sallaa %>%
  ggplot(aes(age.tot, length, color = age.type)) +
  geom_point() +
  facet_wrap(~age.type)

The smolt.only group contains alot of questionable data, i.e. too long fish for being young smolts. Maybe these should be sea age instead

How many are thses fish and where do they come from…

sallaa %>%
  filter(age.type == "smolt.only",
         length > 250) %>%
  ggplot(aes(age.tot, length, color = country)) +
  geom_point() +
  facet_wrap(~age.type)

filter: removed 125,974 rows (>99%), 159 rows remaining

# this is about 160 individuals
sallaa %>%
  filter(age.type == "smolt.only",
         length > 250) %>%
  count()

filter: removed 125,974 rows (>99%), 159 rows remaining
count: now one row and one column, ungrouped

# A tibble: 1 × 1
      n
  <int>
1   159

# if we remove these, the counts per age.type are:
sallaa %>%
  filter(!(age.type == "smolt.only" & length > 250)) %>%
  count(age.type)

filter: removed 159 rows (<1%), 125,974 rows remaining
count: now 3 rows and 2 columns, ungrouped

# A tibble: 3 × 2
  age.type        n
  <chr>       <int>
1 both       106072
2 sea.only    12447
3 smolt.only   7455

# The "smolt.only" comes almost exclusively from Sweden (and some from Finland).
sallaa %>% 
  filter(age.type == "smolt.only") %>%
  count(country)

filter: removed 118,519 rows (94%), 7,614 rows remaining
count: now 2 rows and 2 columns, ungrouped

# A tibble: 2 × 2
  country     n
  <chr>   <int>
1 FIN        62
2 SWE      7552

For now I will use the “both” type only.

# Create a new dataset and reduce number of obs to make models faster
sallaa2 <- sallaa %>%
  filter(age.type == "both") %>% 
  slice_sample(n = 40000)

filter: removed 20,061 rows (16%), 106,072 rows remaining
slice_sample: removed 66,072 rows (62%), 40,000 rows remaining

sallaa2 %>%
  drop_na(sex) %>%
  ggplot(aes(age.tot, length, color = sex)) +
  geom_point() +
  facet_wrap(~sex) +
  scale_x_continuous(breaks = seq(0, 13, 1) )

drop_na: removed 13,544 rows (34%), 26,456 rows remaining

# how old can a salmon get? Is 12 possible or have smolt age been included in sea age?
sallaa2 %>%
  filter(age.tot > 10)

filter: removed 39,995 rows (>99%), 5 rows remaining

# A tibble: 5 × 20
  country  year site       origin length weight age.sea age.sm sex     lat   lon
  <chr>   <dbl> <chr>      <chr>   <dbl>  <dbl>   <dbl>  <dbl> <chr> <dbl> <dbl>
1 FIN      2014 Tornionjo… wild     1100  16500       8      4 f      65.9  24.1
2 FIN      2008 Tornionjo… wild     1190  16600       8      3 f      65.9  24.1
3 SWE      2004 Östersjön… wild     1070 960000       7      4 <NA>   58.5  19.8
4 FIN      2012 Tornionjo… wild     1020   9640       6      5 f      65.9  24.1
5 FIN      2009 Tornionjo… wild     1230  16400       9      3 f      65.9  24.1
# ℹ 9 more variables: asses.unit <dbl>, stock.origin <chr>, stock.unit <chr>,
#   region <chr>, river <chr>, region.gen <chr>, spat.unit <chr>,
#   age.type <chr>, age.tot <dbl>

Hists of age and length

# If dropping all inds where both sea and sm age is NA 
sallaa2 %>% 
  ggplot() +
  geom_density(aes(x = age.tot), fill = "deeppink", alpha = 0.5) +
  xlim(0, 13) +
  
sallaa2 %>% 
  ggplot() + 
  geom_density(aes(x = length), fill = "deeppink", alpha = 0.5)

sallaa2 %>% 
  ggplot() +
  geom_density(aes(x = length, fill = factor(age.tot)), alpha = 0.5) +
  #geom_histogram(aes(x = length)) +
  facet_wrap(~age.tot, scales = "free")

NOTE in code below: Ive silenced warmings to make html manageble when “logProb of data node length[1:nobs]: logProb less than -1e12.”

1. Spatial model w Multivariate normal prior

a. Spatial only

vbgf1a.code<- nimbleCode({
  
  # likelihood
  for(i in 1:nobs){
    length[i] ~ dnorm(mu[i], sd = sig)
    mu[i] <- klt[spat.unit[i],1]*(1-exp(-klt[spat.unit[i],2]*(age.tot[i]-klt[spat.unit[i],3])))
    
    }
  
  # priors 
  sig ~ dgamma(0.01, 0.01)
  
  for (j in 1:nsu){
    klt[j,1:npars] ~ dmnorm(mu_klt[1:npars], tau_klt[1:npars, 1:npars])
  }
  
  mu_klt[1] ~ dnorm(0, 0.01)
  mu_klt[2] ~ dnorm(0, 0.01)
  mu_klt[3] ~ dnorm(0, 0.01)
    
  tau_klt[1:npars, 1:npars] <- inverse(sigma_klt[1:npars, 1:npars])
  
  for(l in 1:npars){
    for(k in 1:npars){
        sigma_klt[l,k] <- Rnew[l,k] * sigma_gpar[l] * sigma_gpar[k]
    }
  }
  
  sigma_gpar[1] ~ dlnorm(0,1)
  sigma_gpar[2] ~ dlnorm(0,1)
  sigma_gpar[3] ~ dlnorm(0,1)

  Rnew[1:npars,1:npars] <- t(R[1:npars,1:npars]) %*% R[1:npars,1:npars]
  alpha[1] <- eta + (npars - 2)/2
  corY[1] ~ dbeta(alpha[1], alpha[1])
  r12 <- 2 * corY[1] - 1

  R[1,1] <- 1
  R[1,2] <- r12
  R[2,2] <- sqrt(1 - r12^2)
  
  #### comment out section below until #### when running with 2 stocks
  R[2:npars,1] <- 0    #with > 2 pars (plus lines below)

  for (m in 2:(npars-1)) {
    ## Draw beta random variable
    alpha[m] <- alpha[(m-1)] - 0.5
    corY[m] ~ dbeta(m / 2, alpha[m])
    ## Draw uniformly on a hypersphere
    for (jj in 1:m) {
      corZ[m, jj] ~ dnorm(0, 1)
    }
    scZ[m, 1:m] <- corZ[m, 1:m] / sqrt(inprod(corZ[m, 1:m], corZ[m, 1:m]))
    R[1:m,(m+1)] <- sqrt(corY[m]) * scZ[m,1:m]
    R[(m+1),(m+1)] <- sqrt(1 - corY[m])
    for(jk in (m+1):npars){
      R[jk,m] <- 0
    }
    
    }  
  })

npars <- 3
inits<-function(){
list(mu_klt = c(rnorm(1,1000,100),
                rnorm(1,0.5,0.5),
                rnorm(1,-1,0.5)),
     sigma_gpar=c(rnorm(1,100,10), rnorm(1,.1,1), rnorm(1,exp(.1),exp(0.05))), corY = c(rbeta(1,2,2),rbeta(1,2,2)))}
# create NIMBLE model
model_3a <- nimbleModel(vbgf1a.code, constants = list(npars=npars, eta=2, nsu = 11,
                                                      nobs = nrow(sallaa2),
                                                      spat.unit = as.integer(factor(sallaa2$spat.unit))),
                        inits=inits(),
                        data = sallaa2 %>% select(age.tot,length) )

select: dropped 18 variables (country, year, site, origin, weight, …)
Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Error in chol.default(model$tau_klt[1:3, 1:3]) : 
  the leading minor of order 1 is not positive


Checking model sizes and dimensions

  [Note] This model is not fully initialized. This is not an error.
         To see which variables are not initialized, use model$initializeInfo().
         For more information on model initialization, see help(modelInitialization).

# compile model
vbgf1a.c <- compileNimble(model_3a)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

# configure mcmc
vbgf1a.conf <- configureMCMC(model_3a, print=TRUE, useConjugacy = FALSE, monitors = "klt")

===== Monitors =====
thin = 1: klt
===== Samplers =====
RW_block sampler (11)
  - klt[]  (11 multivariate elements)
RW sampler (11)
  - sig
  - mu_klt[]  (3 elements)
  - sigma_gpar[]  (3 elements)
  - corZ[]  (2 elements)
  - corY[]  (2 elements)

# build mcmc
vbgf1a.mcmc <- buildMCMC(vbgf1a.conf)
# compile mcmc and specify the project model
vbgf1a.mcmc.c <- compileNimble(vbgf1a.mcmc, project = model_3a)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

# sample mode

Samples

removing output: “warning: logProb of data node length[7]: logProb less than -1e12.”

vbgf1a.samp <- runMCMC(vbgf1a.mcmc.c, nchains = 2, niter = 70000, nburnin = 60000, thin=2, samplesAsCodaMCMC = TRUE)

# NAs produced
summary(vbgf1a.samp)

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 2 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean        SD  Naive SE Time-series SE
klt[1, 1]  1258.23086 3.634e-02 3.634e-04      1.101e-02
klt[2, 1]  1049.07531 1.326e+02 1.326e+00      6.786e+01
klt[3, 1]  1322.23530 5.107e-01 5.107e-03      4.901e-02
klt[4, 1]   964.68721 1.507e+01 1.507e-01      4.951e-01
klt[5, 1]   801.28667 2.012e+01 2.012e-01      5.397e-01
klt[6, 1]   958.28255 1.352e+01 1.352e-01      4.179e-01
klt[7, 1]  1261.03103 1.511e-01 1.511e-03      3.145e-02
klt[8, 1]  1324.38655 2.453e-01 2.453e-03      2.762e-03
klt[9, 1]  1267.53421 3.006e-01 3.006e-03      8.306e-03
klt[10, 1] 1329.55173 2.841e-01 2.841e-03      2.190e-02
klt[11, 1] 1313.42565 5.379e-01 5.379e-03      1.807e-01
klt[1, 2]     0.17771 3.055e-03 3.055e-05      4.980e-04
klt[2, 2]     0.79467 3.837e-01 3.837e-03      1.357e-01
klt[3, 2]     0.29053 2.284e-03 2.284e-05      3.699e-04
klt[4, 2]     0.49711 3.020e-02 3.020e-04      9.555e-04
klt[5, 2]     1.14975 2.545e-01 2.545e-03      1.821e-02
klt[6, 2]     0.85137 1.075e-01 1.075e-03      6.275e-03
klt[7, 2]     0.19670 2.549e-03 2.549e-05      1.013e-03
klt[8, 2]     0.10893 7.063e-03 7.063e-05      4.345e-05
klt[9, 2]     0.16386 1.204e-03 1.204e-05      3.862e-04
klt[10, 2]    0.11673 6.101e-03 6.101e-05      1.281e-03
klt[11, 2]    0.19625 6.732e-03 6.732e-05      1.093e-03
klt[1, 3]    -2.28778 8.706e-02 8.706e-04      1.666e-02
klt[2, 3]     0.57584 9.366e-01 9.366e-03      2.741e-01
klt[3, 3]     1.67222 2.613e-02 2.613e-04      4.686e-03
klt[4, 3]     1.46865 9.160e-02 9.160e-04      2.642e-03
klt[5, 3]     1.28491 5.143e-01 5.143e-03      3.646e-02
klt[6, 3]     1.08142 2.528e-01 2.528e-03      1.733e-02
klt[7, 3]     0.01958 6.326e-02 6.326e-04      2.629e-02
klt[8, 3]    -3.82734 4.222e-01 4.222e-03      3.461e-03
klt[9, 3]    -2.24480 3.283e-02 3.283e-04      1.180e-02
klt[10, 3]   -2.72323 3.363e-01 3.363e-03      7.646e-02
klt[11, 3]   -1.42371 1.238e-01 1.238e-03      2.206e-02

2. Quantiles for each variable:

                 2.5%        25%        50%        75%     97.5%
klt[1, 1]   1.258e+03 1258.20543 1258.24016 1258.25697 1258.2861
klt[2, 1]   9.347e+02  948.75248  959.32516 1212.92276 1264.2492
klt[3, 1]   1.322e+03 1321.69493 1322.31335 1322.72968 1322.8062
klt[4, 1]   9.367e+02  954.23167  964.46026  974.48212  995.2256
klt[5, 1]   7.616e+02  788.11364  801.34608  814.74417  841.7434
klt[6, 1]   9.324e+02  948.94498  958.08843  967.24768  985.9055
klt[7, 1]   1.261e+03 1260.89048 1260.96462 1261.18340 1261.2670
klt[8, 1]   1.324e+03 1324.14230 1324.38504 1324.63253 1324.6428
klt[9, 1]   1.267e+03 1267.24152 1267.52918 1267.82933 1267.8853
klt[10, 1]  1.329e+03 1329.28784 1329.56658 1329.81470 1329.9250
klt[11, 1]  1.313e+03 1312.92511 1313.29751 1313.90467 1314.3019
klt[1, 2]   1.728e-01    0.17459    0.17831    0.18045    0.1822
klt[2, 2]   2.460e-01    0.30427    0.99192    1.10607    1.2665
klt[3, 2]   2.861e-01    0.28895    0.29060    0.29209    0.2950
klt[4, 2]   4.389e-01    0.47607    0.49664    0.51685    0.5575
klt[5, 2]   7.035e-01    0.96518    1.13745    1.30913    1.6875
klt[6, 2]   6.472e-01    0.77973    0.85075    0.92061    1.0694
klt[7, 2]   1.923e-01    0.19442    0.19694    0.19892    0.2007
klt[8, 2]   1.016e-01    0.10187    0.10891    0.11598    0.1164
klt[9, 2]   1.616e-01    0.16294    0.16386    0.16482    0.1660
klt[10, 2]  1.096e-01    0.11084    0.11542    0.12213    0.1276
klt[11, 2]  1.856e-01    0.19112    0.19546    0.20066    0.2112
klt[1, 3]  -2.425e+00   -2.37845   -2.26672   -2.20872   -2.1649
klt[2, 3]  -1.129e+00   -0.50260    1.15041    1.30088    1.4782
klt[3, 3]   1.621e+00    1.65277    1.67362    1.68980    1.7233
klt[4, 3]   1.281e+00    1.40892    1.47237    1.53153    1.6379
klt[5, 3]  -5.307e-03    1.02819    1.38438    1.65064    2.0043
klt[6, 3]   4.863e-01    0.94416    1.11967    1.26044    1.4700
klt[7, 3]  -9.133e-02   -0.03871    0.02786    0.07777    0.1141
klt[8, 3]  -4.266e+00   -4.25032   -3.82819   -3.40645   -3.3841
klt[9, 3]  -2.304e+00   -2.27157   -2.24340   -2.21587   -2.1899
klt[10, 3] -3.102e+00   -3.06305   -2.77596   -2.40935   -2.1524
klt[11, 3] -1.629e+00   -1.51798   -1.43506   -1.33745   -1.1637

mcmc_trace(vbgf1a.samp)

vbgf1a.samp %>%
  gather_draws(klt[spat.unit,par], sep = ",") %>%
  #mutate(val = exp(.value)) %>%
  ggplot() + 
  geom_density(aes(x = .value, color = factor(par))) + # age parameter i
  facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
  theme_light() 

b. Spatial w sex covar. for Linf

vbgf1b.code <- nimbleCode({
  
  # likelihood
  for(i in 1:nobs){
    length[i] ~ dnorm(mu[i], sd = sig)
    mu[i] <- (klt[spat.unit[i],1] + beta1*sex[i]) * (1-exp(-klt[spat.unit[i],2]*(age.tot[i]-klt[spat.unit[i],3])))
    #mu[i] <- (klt[spat.unit[i],1] + beta1*sex[i]) * (1-exp(-(klt[spat.unit[i],2] + beta1*sex[i])*(age.tot[i]-(klt[spat.unit[i],3] + beta1*sex[i]))))
    }
  
  # priors 
  sig ~ dgamma(0.01, 0.01)
  
  beta1 ~ dnorm(0, 0.01)
  
  for (j in 1:nsu){
    klt[j,1:npars] ~ dmnorm(mu_klt[1:npars], tau_klt[1:npars, 1:npars])
  }
  
  mu_klt[1] ~ dnorm(0, 0.01)
  mu_klt[2] ~ dnorm(0, 0.01)
  mu_klt[3] ~ dnorm(0, 0.01)
    
  tau_klt[1:npars, 1:npars] <- inverse(sigma_klt[1:npars, 1:npars])
  
  for(l in 1:npars){
    for(k in 1:npars){
        sigma_klt[l,k] <- Rnew[l,k] * sigma_gpar[l] * sigma_gpar[k]
    }
  }
  
  sigma_gpar[1] ~ dlnorm(0,1)
  sigma_gpar[2] ~ dlnorm(0,1)
  sigma_gpar[3] ~ dlnorm(0,1)

  Rnew[1:npars,1:npars] <- t(R[1:npars,1:npars]) %*% R[1:npars,1:npars]
  alpha[1] <- eta + (npars - 2)/2
  corY[1] ~ dbeta(alpha[1], alpha[1])
  r12 <- 2 * corY[1] - 1

  R[1,1] <- 1
  R[1,2] <- r12
  R[2,2] <- sqrt(1 - r12^2)
  
  #### comment out section below until #### when running with 2 stocks
  R[2:npars,1] <- 0    #with > 2 pars (plus lines below)

  for (m in 2:(npars-1)) {
    ## Draw beta random variable
    alpha[m] <- alpha[(m-1)] - 0.5
    corY[m] ~ dbeta(m / 2, alpha[m])
    ## Draw uniformly on a hypersphere
    for (jj in 1:m) {
      corZ[m, jj] ~ dnorm(0, 1)
    }
    scZ[m, 1:m] <- corZ[m, 1:m] / sqrt(inprod(corZ[m, 1:m], corZ[m, 1:m]))
    R[1:m,(m+1)] <- sqrt(corY[m]) * scZ[m,1:m]
    R[(m+1),(m+1)] <- sqrt(1 - corY[m])
    for(jk in (m+1):npars){
      R[jk,m] <- 0
    }
    
    }  
  })

npars <- 3
inits<-function(){
list(mu_klt = c(rnorm(1,1100,100),
                rnorm(1,0.5,0.5),
                rnorm(1,-1,0.5)),
     beta1 = rnorm(1,-15,5),
     sigma_gpar=c(rnorm(1,100,10), rnorm(1,.1,1), rnorm(1,exp(.1),exp(0.05))), corY = c(rbeta(1,2,2),rbeta(1,2,2)))}

data_3b = sallaa2 %>% select(age.tot, length, sex, spat.unit) %>% 
  mutate(sex = as.integer(if_else(sex == "f", 0, 1)),
         spat.unit = as.integer(factor(sallaa2$spat.unit))) %>% drop_na(sex)

select: dropped 16 variables (country, year, site, origin, weight, …)
mutate: converted 'sex' from character to integer (0 new NA)
        converted 'spat.unit' from character to integer (0 new NA)
drop_na: removed 13,544 rows (34%), 26,456 rows remaining

# create NIMBLE model
model_3b <- nimbleModel(vbgf1b.code, constants = list(npars=npars, eta=2, nsu = 11,
                                                      nobs = nrow(data_3b),
                                                      spat.unit = data_3b$spat.unit),
                        inits=inits(),
                        data = data_3b %>% select(age.tot,length,sex))

select: dropped one variable (spat.unit)
Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Error in chol.default(model$tau_klt[1:3, 1:3]) : 
  the leading minor of order 1 is not positive


Checking model sizes and dimensions

  [Note] This model is not fully initialized. This is not an error.
         To see which variables are not initialized, use model$initializeInfo().
         For more information on model initialization, see help(modelInitialization).

# compile model
vbgf1b.c <- compileNimble(model_3b)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

vbgf1b.c$beta1 # if not using inits in nimbleModel, these calculated varaibles wont show up

[1] -16.90863

vbgf1b.c$mu_klt

[1] 870.7926150   0.5817883  -1.3370027

# configure mcmc
vbgf1b.conf <- configureMCMC(vbgf1b.c, print=TRUE, useConjugacy = FALSE, , monitors =  c("klt", "beta1"))

===== Monitors =====
thin = 1: beta1, klt
===== Samplers =====
RW_block sampler (11)
  - klt[]  (11 multivariate elements)
RW sampler (12)
  - sig
  - beta1
  - mu_klt[]  (3 elements)
  - sigma_gpar[]  (3 elements)
  - corZ[]  (2 elements)
  - corY[]  (2 elements)

# build mcmc
vbgf1b.mcmc <- buildMCMC(vbgf1b.conf)
# compile mcmc and specify the project model
vbgf1b.mcmc.c <- compileNimble(vbgf1b.mcmc)#, project = vbgf1b.c)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

# sample mode
vbgf1b.conf$getUnsampledNodes()

character(0)

Samples

removing output: “warning: logProb of data node length[7]: logProb less than -1e12.”

vbgf1b.samp <- runMCMC(vbgf1b.mcmc.c, nchains = 2, niter = 150000, nburnin = 140000, thin=2, samplesAsCodaMCMC = TRUE)

summary(vbgf1b.samp)

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 2 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean       SD  Naive SE Time-series SE
beta1       -13.4757  1.32568 0.0132568      0.0334847
klt[1, 1]   847.8016  7.70717 0.0770717      0.2499661
klt[2, 1]   912.0386  7.96798 0.0796798      0.1973302
klt[3, 1]  1115.4382  6.10837 0.0610837      0.1651783
klt[4, 1]   894.0155 10.51094 0.1051094      0.2566106
klt[5, 1]   717.5660 16.25467 0.1625467      0.5200636
klt[6, 1]   895.6693 10.22345 0.1022345      0.2593153
klt[7, 1]   986.5423  8.94990 0.0894990      0.2302578
klt[8, 1]   803.0491  8.82576 0.0882576      0.6906068
klt[9, 1]   795.2545  8.72346 0.0872346      0.3886411
klt[10, 1]  843.8859  9.19154 0.0919154      0.2370282
klt[11, 1]  678.3877  2.66195 0.0266195      0.0815615
klt[1, 2]     0.9267  0.07345 0.0007345      0.0024529
klt[2, 2]     1.7037  0.19801 0.0019801      0.0051776
klt[3, 2]     0.5080  0.01134 0.0001134      0.0003017
klt[4, 2]     0.6602  0.03823 0.0003823      0.0009821
klt[5, 2]     1.6197  0.59472 0.0059472      0.0216279
klt[6, 2]     1.3668  0.20428 0.0020428      0.0060324
klt[7, 2]     0.4807  0.02031 0.0002031      0.0005243
klt[8, 2]     0.5944  0.04455 0.0004455      0.0034444
klt[9, 2]     1.0707  0.12458 0.0012458      0.0064316
klt[10, 2]    0.5324  0.03460 0.0003460      0.0008995
klt[11, 2]   -3.2863  1.11230 0.0111230      0.0592674
klt[1, 3]     0.3194  0.10502 0.0010502      0.0034707
klt[2, 3]     1.4032  0.14342 0.0014342      0.0036059
klt[3, 3]     2.2692  0.03064 0.0003064      0.0008041
klt[4, 3]     1.6555  0.08498 0.0008498      0.0020868
klt[5, 3]    -0.0184  1.16712 0.0116712      0.0436991
klt[6, 3]     1.2311  0.25346 0.0025346      0.0068775
klt[7, 3]     1.4998  0.07880 0.0007880      0.0019685
klt[8, 3]    -0.5953  0.13634 0.0013634      0.0094398
klt[9, 3]     0.4233  0.14048 0.0014048      0.0069024
klt[10, 3]    0.1218  0.14015 0.0014015      0.0035288
klt[11, 3]    6.8230  0.88065 0.0088065      0.0344810

2. Quantiles for each variable:

                2.5%        25%       50%       75%     97.5%
beta1       -16.0140  -14.36563  -13.4727  -12.5911  -10.8911
klt[1, 1]   833.5547  842.60992  847.3445  852.8821  863.5855
klt[2, 1]   896.9981  906.48934  911.8794  917.2570  928.3083
klt[3, 1]  1103.5217 1111.21230 1115.4199 1119.4835 1127.6751
klt[4, 1]   873.2189  886.99645  894.0081  900.9939  915.1755
klt[5, 1]   686.2695  706.50716  717.7139  728.2246  750.1994
klt[6, 1]   876.3707  888.79613  895.3514  902.3616  917.0833
klt[7, 1]   969.3778  980.24544  986.5072  992.4844 1004.7412
klt[8, 1]   786.8205  796.64950  802.6614  808.7574  821.5133
klt[9, 1]   778.7407  789.33994  794.9916  801.1061  813.5274
klt[10, 1]  826.9499  837.59273  843.6372  849.8386  862.9392
klt[11, 1]  673.2730  676.58769  678.3695  680.1730  683.6762
klt[1, 2]     0.7920    0.87656    0.9240    0.9745    1.0755
klt[2, 2]     1.3570    1.56500    1.6910    1.8259    2.1349
klt[3, 2]     0.4860    0.50023    0.5078    0.5156    0.5304
klt[4, 2]     0.5898    0.63415    0.6584    0.6852    0.7388
klt[5, 2]     0.8248    1.19985    1.4895    1.9039    3.1254
klt[6, 2]     1.0237    1.22800    1.3475    1.4843    1.8318
klt[7, 2]     0.4424    0.46721    0.4800    0.4942    0.5215
klt[8, 2]     0.5084    0.56375    0.5942    0.6241    0.6846
klt[9, 2]     0.8612    0.98303    1.0600    1.1434    1.3690
klt[10, 2]    0.4655    0.50897    0.5320    0.5555    0.6044
klt[11, 2]   -5.9285   -3.94054   -3.1291   -2.5098   -1.5762
klt[1, 3]     0.1033    0.25047    0.3209    0.3927    0.5130
klt[2, 3]     1.0676    1.31644    1.4213    1.5061    1.6320
klt[3, 3]     2.2089    2.24817    2.2688    2.2903    2.3292
klt[4, 3]     1.4834    1.60084    1.6588    1.7135    1.8145
klt[5, 3]    -2.6866   -0.72943    0.1377    0.8196    1.8712
klt[6, 3]     0.6188    1.09212    1.2731    1.4113    1.6051
klt[7, 3]     1.3433    1.44896    1.5005    1.5533    1.6511
klt[8, 3]    -0.8828   -0.68550   -0.5860   -0.4992   -0.3462
klt[9, 3]     0.1467    0.33063    0.4245    0.5148    0.7129
klt[10, 3]   -0.1745    0.03592    0.1264    0.2184    0.3796
klt[11, 3]    5.6773    6.18279    6.6224    7.2599    9.0017

mcmc_trace(vbgf1b.samp)

gelman.diag(vbgf1b.samp)

Potential scale reduction factors:

           Point est. Upper C.I.
beta1            1.00       1.00
klt[1, 1]        1.00       1.01
klt[2, 1]        1.00       1.00
klt[3, 1]        1.00       1.00
klt[4, 1]        1.01       1.02
klt[5, 1]        1.00       1.00
klt[6, 1]        1.00       1.00
klt[7, 1]        1.01       1.05
klt[8, 1]        1.04       1.16
klt[9, 1]        1.00       1.00
klt[10, 1]       1.01       1.04
klt[11, 1]       1.00       1.00
klt[1, 2]        1.00       1.01
klt[2, 2]        1.00       1.00
klt[3, 2]        1.00       1.02
klt[4, 2]        1.00       1.00
klt[5, 2]        1.00       1.00
klt[6, 2]        1.00       1.01
klt[7, 2]        1.00       1.02
klt[8, 2]        1.04       1.16
klt[9, 2]        1.01       1.01
klt[10, 2]       1.02       1.08
klt[11, 2]       1.18       1.62
klt[1, 3]        1.00       1.02
klt[2, 3]        1.00       1.00
klt[3, 3]        1.01       1.03
klt[4, 3]        1.00       1.00
klt[5, 3]        1.04       1.14
klt[6, 3]        1.00       1.01
klt[7, 3]        1.00       1.00
klt[8, 3]        1.03       1.13
klt[9, 3]        1.00       1.00
klt[10, 3]       1.02       1.08
klt[11, 3]       1.41       2.50

Multivariate psrf

1.25

effectiveSize(vbgf1b.samp)

     beta1  klt[1, 1]  klt[2, 1]  klt[3, 1]  klt[4, 1]  klt[5, 1]  klt[6, 1] 
 1587.7599   959.9678  1630.3322  1370.0104  1691.4868  1022.1972  1563.3609 
 klt[7, 1]  klt[8, 1]  klt[9, 1] klt[10, 1] klt[11, 1]  klt[1, 2]  klt[2, 2] 
 1521.3350   388.3966   519.3222  1516.0377  1071.5903   896.1996  1466.4085 
 klt[3, 2]  klt[4, 2]  klt[5, 2]  klt[6, 2]  klt[7, 2]  klt[8, 2]  klt[9, 2] 
 1414.0046  1545.8230   790.7120  1174.9287  1512.4240   387.3523   456.7571 
klt[10, 2] klt[11, 2]  klt[1, 3]  klt[2, 3]  klt[3, 3]  klt[4, 3]  klt[5, 3] 
 1491.3645   372.7850   914.8764  1589.8267  1452.7755  1665.3688   699.4474 
 klt[6, 3]  klt[7, 3]  klt[8, 3]  klt[9, 3] klt[10, 3] klt[11, 3] 
 1384.1219  1626.5421   431.7080   499.6589  1573.2031   622.4520 

vbgf1b.samp %>%
  gather_draws(klt[spat.unit,par], beta1, sep = ",") %>%
  ggplot() + 
  geom_density(aes(x = .value, color = factor(par))) + # age parameter i
  facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
  theme_light() 

c. Spatial w sex covar. for all vbgf pars

vbgf1c.code <- nimbleCode({
  
  # likelihood
  for(i in 1:nobs){
    length[i] ~ dnorm(mu[i], sd = sig)
    mu[i] <- (klt[spat.unit[i],1] + beta1*sex[i]) * (1-exp(-(klt[spat.unit[i],2] + beta2*sex[i])*(age.tot[i]-(klt[spat.unit[i],3] + beta3*sex[i]))))
    }
  
  # priors 
  sig ~ dgamma(0.01, 0.01)
  
  beta1 ~ dnorm(0, 0.01)
  beta2 ~ dnorm(0, 0.01)
  beta3 ~ dnorm(0, 0.01)
  
  for (j in 1:nsu){
    klt[j,1:npars] ~ dmnorm(mu_klt[1:npars], tau_klt[1:npars, 1:npars])
  }
  
  mu_klt[1] ~ dnorm(0, 0.01)
  mu_klt[2] ~ dnorm(0, 0.01)
  mu_klt[3] ~ dnorm(0, 0.01)
    
  tau_klt[1:npars, 1:npars] <- inverse(sigma_klt[1:npars, 1:npars])
  
  for(l in 1:npars){
    for(k in 1:npars){
        sigma_klt[l,k] <- Rnew[l,k] * sigma_gpar[l] * sigma_gpar[k]
    }
  }
  
  sigma_gpar[1] ~ dlnorm(0,1)
  sigma_gpar[2] ~ dlnorm(0,1)
  sigma_gpar[3] ~ dlnorm(0,1)

  Rnew[1:npars,1:npars] <- t(R[1:npars,1:npars]) %*% R[1:npars,1:npars]
  alpha[1] <- eta + (npars - 2)/2
  corY[1] ~ dbeta(alpha[1], alpha[1])
  r12 <- 2 * corY[1] - 1

  R[1,1] <- 1
  R[1,2] <- r12
  R[2,2] <- sqrt(1 - r12^2)
  
  #### comment out section below until #### when running with 2 stocks
  R[2:npars,1] <- 0    #with > 2 pars (plus lines below)

  for (m in 2:(npars-1)) {
    ## Draw beta random variable
    alpha[m] <- alpha[(m-1)] - 0.5
    corY[m] ~ dbeta(m / 2, alpha[m])
    ## Draw uniformly on a hypersphere
    for (jj in 1:m) {
      corZ[m, jj] ~ dnorm(0, 1)
    }
    scZ[m, 1:m] <- corZ[m, 1:m] / sqrt(inprod(corZ[m, 1:m], corZ[m, 1:m]))
    R[1:m,(m+1)] <- sqrt(corY[m]) * scZ[m,1:m]
    R[(m+1),(m+1)] <- sqrt(1 - corY[m])
    for(jk in (m+1):npars){
      R[jk,m] <- 0
    }
    
    }  
  })

npars <- 3
inits<-function(){
list(mu_klt = c(rnorm(1,1000,100),
                rnorm(1,0.5,0.5),
                rnorm(1,0,1)),
     beta1 = rnorm(1,100,50),
     beta2 = rnorm(1,0.5,0.5),
     beta3 = rnorm(1,0,1),
     sigma_gpar=c(rnorm(1,100,10), rnorm(1,.1,1), rnorm(1,exp(.1),exp(0.05))), corY = c(rbeta(1,2,2),rbeta(1,2,2)))}

data_3d = sallaa2 %>% select(age.tot, length, sex, spat.unit) %>% 
  mutate(sex = as.integer(if_else(sex == "f", 0, 1)),
         spat.unit = as.integer(factor(sallaa2$spat.unit))) %>% drop_na(sex)

select: dropped 16 variables (country, year, site, origin, weight, …)
mutate: converted 'sex' from character to integer (0 new NA)
        converted 'spat.unit' from character to integer (0 new NA)
drop_na: removed 13,544 rows (34%), 26,456 rows remaining

# create NIMBLE model
model_3d <- nimbleModel(vbgf1c.code, constants = list(npars=npars, eta=2, nsu = 11,
                                                      nobs = nrow(data_3b),
                                                      spat.unit = data_3b$spat.unit),
                        inits=inits(),
                        data = data_3b %>% select(age.tot,length,sex))

select: dropped one variable (spat.unit)
Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Error in chol.default(model$tau_klt[1:3, 1:3]) : 
  the leading minor of order 1 is not positive


Checking model sizes and dimensions

  [Note] This model is not fully initialized. This is not an error.
         To see which variables are not initialized, use model$initializeInfo().
         For more information on model initialization, see help(modelInitialization).

# compile model
vbgf1c.c <- compileNimble(model_3d)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

# configure mcmc
vbgf1c.conf <- configureMCMC(vbgf1c.c, print=TRUE, useConjugacy = FALSE, , monitors =  c("klt", "beta1", "beta2", "beta3"), enableWAIC = TRUE)

===== Monitors =====
thin = 1: beta1, beta2, beta3, klt
===== Samplers =====
RW_block sampler (11)
  - klt[]  (11 multivariate elements)
RW sampler (14)
  - sig
  - beta1
  - beta2
  - beta3
  - mu_klt[]  (3 elements)
  - sigma_gpar[]  (3 elements)
  - corZ[]  (2 elements)
  - corY[]  (2 elements)

# build mcmc
vbgf1c.mcmc <- buildMCMC(vbgf1c.conf)
# compile mcmc and specify the project model
vbgf1c.mcmc.c <- compileNimble(vbgf1c.mcmc)#, project = vbgf1b.c)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Samples

removing output: “warning: logProb of data node length[7]: logProb less than -1e12.”

# sample model
vbgf1c.samp <- runMCMC(vbgf1c.mcmc.c, nchains = 2, niter = 250000, nburnin = 240000, thin=2, samplesAsCodaMCMC = TRUE, WAIC=TRUE)
# vbgf1c.samp <- runMCMC(vbgf1c.mcmc.c, nchains = 2, niter = 25000, nburnin = 24000, thin=2, samplesAsCodaMCMC = TRUE, WAIC=TRUE)

summary(vbgf1c.samp$samples)

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 2 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean        SD  Naive SE Time-series SE
beta1       9.463e+01 6.188e+00 6.188e-02      5.830e-01
beta2      -6.381e-03 1.359e-03 1.359e-05      3.257e-04
beta3       6.076e-01 4.289e-02 4.289e-04      4.946e-03
klt[1, 1]   1.060e+03 2.897e+01 2.897e-01      1.838e+00
klt[2, 1]   1.105e+03 5.230e+01 5.230e-01      4.081e+00
klt[3, 1]   1.192e+03 1.057e+01 1.057e-01      3.543e-01
klt[4, 1]   2.401e+03 2.059e+02 2.059e+00      4.140e+01
klt[5, 1]   1.348e+03 2.880e+02 2.880e+00      4.565e+01
klt[6, 1]   1.186e+03 8.067e+01 8.067e-01      5.398e+00
klt[7, 1]   1.790e+03 1.030e+02 1.030e+00      2.164e+01
klt[8, 1]   1.195e+03 4.336e+01 4.336e-01      1.241e+01
klt[9, 1]   1.069e+03 3.450e+01 3.450e-01      2.590e+00
klt[10, 1]  1.067e+03 4.326e+01 4.326e-01      2.739e+00
klt[11, 1]  9.184e+02 3.201e+01 3.201e-01      1.918e+00
klt[1, 2]   2.559e-01 2.441e-02 2.441e-04      1.627e-03
klt[2, 2]   4.574e-01 1.182e-01 1.182e-03      1.154e-02
klt[3, 2]   3.462e-01 1.043e-02 1.043e-04      4.082e-04
klt[4, 2]   6.418e-02 9.162e-03 9.162e-05      1.810e-03
klt[5, 2]   2.191e-01 1.029e-01 1.029e-03      1.039e-02
klt[6, 2]   2.964e-01 7.163e-02 7.163e-04      4.467e-03
klt[7, 2]   9.384e-02 1.013e-02 1.013e-04      1.872e-03
klt[8, 2]   1.235e-01 1.072e-02 1.072e-04      2.444e-03
klt[9, 2]   2.129e-01 2.210e-02 2.210e-04      1.516e-03
klt[10, 2]  1.844e-01 2.332e-02 2.332e-04      1.540e-03
klt[11, 2]  3.731e-01 5.071e-02 5.071e-04      3.434e-03
klt[1, 3]  -2.020e+00 2.070e-01 2.070e-03      1.360e-02
klt[2, 3]  -1.409e-01 5.043e-01 5.043e-03      4.381e-02
klt[3, 3]   1.441e+00 6.046e-02 6.046e-04      4.597e-03
klt[4, 3]  -1.855e+00 3.065e-01 3.065e-03      2.920e-02
klt[5, 3]  -1.088e+00 9.969e-01 9.969e-03      1.268e-01
klt[6, 3]  -7.653e-01 5.917e-01 5.917e-03      3.524e-02
klt[7, 3]  -1.507e+00 2.173e-01 2.173e-03      2.260e-02
klt[8, 3]  -4.217e+00 2.389e-01 2.389e-03      2.951e-02
klt[9, 3]  -2.477e+00 2.445e-01 2.445e-03      1.748e-02
klt[10, 3] -2.222e+00 3.590e-01 3.590e-03      2.421e-02
klt[11, 3] -1.398e+00 2.577e-01 2.577e-03      1.783e-02

2. Quantiles for each variable:

                 2.5%        25%        50%        75%      97.5%
beta1       8.173e+01  9.046e+01  9.495e+01  9.889e+01  1.060e+02
beta2      -9.933e-03 -7.025e-03 -6.285e-03 -5.441e-03 -4.183e-03
beta3       5.159e-01  5.817e-01  6.074e-01  6.367e-01  6.879e-01
klt[1, 1]   1.004e+03  1.040e+03  1.060e+03  1.081e+03  1.115e+03
klt[2, 1]   1.009e+03  1.068e+03  1.103e+03  1.139e+03  1.220e+03
klt[3, 1]   1.172e+03  1.185e+03  1.192e+03  1.199e+03  1.213e+03
klt[4, 1]   1.987e+03  2.265e+03  2.384e+03  2.550e+03  2.808e+03
klt[5, 1]   9.913e+02  1.126e+03  1.260e+03  1.509e+03  1.987e+03
klt[6, 1]   1.060e+03  1.126e+03  1.176e+03  1.235e+03  1.375e+03
klt[7, 1]   1.592e+03  1.714e+03  1.785e+03  1.865e+03  1.989e+03
klt[8, 1]   1.104e+03  1.169e+03  1.187e+03  1.222e+03  1.288e+03
klt[9, 1]   1.001e+03  1.044e+03  1.069e+03  1.093e+03  1.134e+03
klt[10, 1]  9.908e+02  1.036e+03  1.065e+03  1.098e+03  1.155e+03
klt[11, 1]  8.572e+02  8.955e+02  9.176e+02  9.401e+02  9.840e+02
klt[1, 2]   2.166e-01  2.374e-01  2.537e-01  2.708e-01  3.108e-01
klt[2, 2]   2.791e-01  3.750e-01  4.394e-01  5.194e-01  7.492e-01
klt[3, 2]   3.257e-01  3.392e-01  3.463e-01  3.532e-01  3.669e-01
klt[4, 2]   4.915e-02  5.727e-02  6.378e-02  6.945e-02  8.565e-02
klt[5, 2]   7.725e-02  1.346e-01  2.053e-01  2.868e-01  4.556e-01
klt[6, 2]   1.786e-01  2.422e-01  2.914e-01  3.433e-01  4.492e-01
klt[7, 2]   7.771e-02  8.587e-02  9.323e-02  1.010e-01  1.150e-01
klt[8, 2]   1.018e-01  1.161e-01  1.243e-01  1.301e-01  1.477e-01
klt[9, 2]   1.759e-01  1.968e-01  2.109e-01  2.271e-01  2.607e-01
klt[10, 2]  1.444e-01  1.672e-01  1.832e-01  2.001e-01  2.316e-01
klt[11, 2]  2.860e-01  3.375e-01  3.686e-01  4.035e-01  4.901e-01
klt[1, 3]  -2.386e+00 -2.173e+00 -2.024e+00 -1.883e+00 -1.585e+00
klt[2, 3]  -1.197e+00 -4.623e-01 -1.122e-01  2.089e-01  7.851e-01
klt[3, 3]   1.326e+00  1.402e+00  1.440e+00  1.480e+00  1.568e+00
klt[4, 3]  -2.475e+00 -2.050e+00 -1.842e+00 -1.645e+00 -1.280e+00
klt[5, 3]  -3.289e+00 -1.791e+00 -9.760e-01 -2.956e-01  4.885e-01
klt[6, 3]  -2.029e+00 -1.159e+00 -7.001e-01 -3.313e-01  2.033e-01
klt[7, 3]  -1.913e+00 -1.670e+00 -1.512e+00 -1.348e+00 -1.099e+00
klt[8, 3]  -4.748e+00 -4.363e+00 -4.198e+00 -4.054e+00 -3.786e+00
klt[9, 3]  -2.940e+00 -2.651e+00 -2.471e+00 -2.306e+00 -2.010e+00
klt[10, 3] -2.939e+00 -2.471e+00 -2.203e+00 -1.963e+00 -1.569e+00
klt[11, 3] -1.944e+00 -1.561e+00 -1.397e+00 -1.227e+00 -8.910e-01

mcmc_trace(vbgf1c.samp$samples)

gelman.diag(vbgf1c.samp$samples)

Potential scale reduction factors:

           Point est. Upper C.I.
beta1            1.09       1.35
beta2            1.25       2.05
beta3            1.03       1.06
klt[1, 1]        1.00       1.00
klt[2, 1]        1.02       1.08
klt[3, 1]        1.00       1.01
klt[4, 1]        1.16       1.57
klt[5, 1]        1.28       2.20
klt[6, 1]        1.01       1.03
klt[7, 1]        1.11       1.28
klt[8, 1]        1.22       1.22
klt[9, 1]        1.00       1.00
klt[10, 1]       1.03       1.14
klt[11, 1]       1.01       1.06
klt[1, 2]        1.00       1.00
klt[2, 2]        1.09       1.21
klt[3, 2]        1.01       1.03
klt[4, 2]        1.21       1.76
klt[5, 2]        1.07       1.21
klt[6, 2]        1.03       1.08
klt[7, 2]        1.13       1.37
klt[8, 2]        1.21       1.25
klt[9, 2]        1.01       1.01
klt[10, 2]       1.04       1.12
klt[11, 2]       1.04       1.13
klt[1, 3]        1.00       1.00
klt[2, 3]        1.01       1.02
klt[3, 3]        1.02       1.05
klt[4, 3]        1.07       1.09
klt[5, 3]        1.12       1.40
klt[6, 3]        1.01       1.03
klt[7, 3]        1.02       1.04
klt[8, 3]        1.12       1.14
klt[9, 3]        1.02       1.02
klt[10, 3]       1.02       1.03
klt[11, 3]       1.03       1.11

Multivariate psrf

1.67

effectiveSize(vbgf1c.samp$samples)

     beta1      beta2      beta3  klt[1, 1]  klt[2, 1]  klt[3, 1]  klt[4, 1] 
 109.54888   56.70203   73.59051  246.40470  186.53019  888.25937   48.66263 
 klt[5, 1]  klt[6, 1]  klt[7, 1]  klt[8, 1]  klt[9, 1] klt[10, 1] klt[11, 1] 
  42.58177  277.45855   28.60732   14.20327  175.01093  293.09047  310.87828 
 klt[1, 2]  klt[2, 2]  klt[3, 2]  klt[4, 2]  klt[5, 2]  klt[6, 2]  klt[7, 2] 
 225.77238  145.96322  660.21354   51.18773  106.85682  287.40404   43.54309 
 klt[8, 2]  klt[9, 2] klt[10, 2] klt[11, 2]  klt[1, 3]  klt[2, 3]  klt[3, 3] 
  27.66751  210.99769  289.51123  268.62571  230.25583  151.97552  182.53055 
 klt[4, 3]  klt[5, 3]  klt[6, 3]  klt[7, 3]  klt[8, 3]  klt[9, 3] klt[10, 3] 
 254.06818   98.85939  316.36529   99.82188  133.07078  192.00393  254.08059 
klt[11, 3] 
 258.30060 

vbgf1c.samp$samples %>%
  gather_draws(klt[spat.unit,par], beta1, sep = ",") %>%
  filter(par == 1) %>%
  ggplot() + 
  geom_density(aes(x = .value, color = factor(par))) + # age parameter i
  facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
  theme_light() 

filter (grouped): removed 230,000 rows (68%), 110,000 rows remaining (removed
23 groups, 11 groups remaining)

vbgf1c.samp$samples %>%
  gather_draws(klt[spat.unit,par], beta1, sep = ",") %>%
  filter(par == 2) %>%
  ggplot() + 
  geom_density(aes(x = .value, color = factor(par))) + # age parameter i
  facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
  theme_light() 

filter (grouped): removed 230,000 rows (68%), 110,000 rows remaining (removed
23 groups, 11 groups remaining)

vbgf1c.samp$samples %>%
  gather_draws(klt[spat.unit,par], beta1, sep = ",") %>%
  filter(par == 3) %>%
  ggplot() + 
  geom_density(aes(x = .value, color = factor(par))) + # age parameter i
  facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
  theme_light() 

filter (grouped): removed 230,000 rows (68%), 110,000 rows remaining (removed
23 groups, 11 groups remaining)

d. Spatial sex covar. for all vbgf pars

vbgf1d.code <- nimbleCode({
  
  # likelihood
  for(i in 1:nobs){
    length[i] ~ dnorm(mu[i], sd = sig)
    mu[i] <-  (klt[spat.unit[i],1] + b1[spat.unit[i]]*sex[i]) *
      (1-exp(-(klt[spat.unit[i],2] + b2[spat.unit[i]]*sex[i]) * (age.tot[i]-(klt[spat.unit[i],3] + b3[spat.unit[i]]*sex[i]))))
    }
  
  # priors 
  sig ~ dgamma(0.01, 0.01)
  
  for (j in 1:nsu){
    klt[j,1:npars] ~ dmnorm(mu_klt[1:npars], tau_klt[1:npars, 1:npars])
    
    b1[j] ~ dnorm(0, 0.01)
    b2[j] ~ dnorm(0, 0.01)
    b3[j] ~ dnorm(0, 0.01)
  }
  
  mu_klt[1] ~ dnorm(0, 0.01)
  mu_klt[2] ~ dnorm(0, 0.01)
  mu_klt[3] ~ dnorm(0, 0.01)
    
  tau_klt[1:npars, 1:npars] <- inverse(sigma_klt[1:npars, 1:npars])
  
  for(l in 1:npars){
    for(k in 1:npars){
        sigma_klt[l,k] <- Rnew[l,k] * sigma_gpar[l] * sigma_gpar[k]
    }
  }
  
  sigma_gpar[1] ~ dlnorm(0,1)
  sigma_gpar[2] ~ dlnorm(0,1)
  sigma_gpar[3] ~ dlnorm(0,1)

  Rnew[1:npars,1:npars] <- t(R[1:npars,1:npars]) %*% R[1:npars,1:npars]
  alpha[1] <- eta + (npars - 2)/2
  corY[1] ~ dbeta(alpha[1], alpha[1])
  r12 <- 2 * corY[1] - 1

  R[1,1] <- 1
  R[1,2] <- r12
  R[2,2] <- sqrt(1 - r12^2)
  
  #### comment out section below until #### when running with 2 stocks
  R[2:npars,1] <- 0    #with > 2 pars (plus lines below)

  for (m in 2:(npars-1)) {
    ## Draw beta random variable
    alpha[m] <- alpha[(m-1)] - 0.5
    corY[m] ~ dbeta(m / 2, alpha[m])
    ## Draw uniformly on a hypersphere
    for (jj in 1:m) {
      corZ[m, jj] ~ dnorm(0, 1)
    }
    scZ[m, 1:m] <- corZ[m, 1:m] / sqrt(inprod(corZ[m, 1:m], corZ[m, 1:m]))
    R[1:m,(m+1)] <- sqrt(corY[m]) * scZ[m,1:m]
    R[(m+1),(m+1)] <- sqrt(1 - corY[m])
    for(jk in (m+1):npars){
      R[jk,m] <- 0
    }
    
    }  
  })

npars <- 3
nsu <- length(unique(sallaa2$spat.unit))

inits<-function(){
list(mu_klt = c(rnorm(1,1000,100),
                rnorm(1,0.5,0.5),
                rnorm(1,0,1)),
     b1 = rnorm(nsu,100,50),
     b2 = rnorm(nsu,0.5,0.5),
     b3 = rnorm(nsu,0,1),
     sigma_gpar=c(rlnorm(1,-0.60,0.30),rlnorm(1,-3,0.05),rlnorm(1,-1,0.05)), corY = c(rbeta(1,2,2),rbeta(1,2,2)))}

data_3d = sallaa2 %>% select(age.tot, length, sex, spat.unit) %>% 
  mutate(sex = as.integer(if_else(sex == "f", 0, 1)),
         spat.unit = as.integer(factor(sallaa2$spat.unit))) %>% drop_na(sex)

select: dropped 16 variables (country, year, site, origin, weight, …)
mutate: converted 'sex' from character to integer (0 new NA)
        converted 'spat.unit' from character to integer (0 new NA)
drop_na: removed 13,544 rows (34%), 26,456 rows remaining

# create NIMBLE model
model_3d <- nimbleModel(vbgf1d.code, constants = list(npars=npars, eta=2, nsu = nsu,
                                                      nobs = nrow(data_3b),
                                                      spat.unit = data_3b$spat.unit),
                        inits=inits(),
                        data = data_3b %>% select(age.tot,length,sex))

select: dropped one variable (spat.unit)
Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Error in chol.default(model$tau_klt[1:3, 1:3]) : 
  the leading minor of order 1 is not positive


Checking model sizes and dimensions

  [Note] This model is not fully initialized. This is not an error.
         To see which variables are not initialized, use model$initializeInfo().
         For more information on model initialization, see help(modelInitialization).

# compile model
vbgf1d.c <- compileNimble(model_3d)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

# configure mcmc
vbgf1d.conf <- configureMCMC(vbgf1d.c, print=TRUE, useConjugacy = FALSE, , monitors =  c("klt", "b1", "b2", "b3"), enableWAIC = TRUE)

===== Monitors =====
thin = 1: b1, b2, b3, klt
===== Samplers =====
RW_block sampler (11)
  - klt[]  (11 multivariate elements)
RW sampler (44)
  - sig
  - b1[]  (11 elements)
  - b2[]  (11 elements)
  - b3[]  (11 elements)
  - mu_klt[]  (3 elements)
  - sigma_gpar[]  (3 elements)
  - corZ[]  (2 elements)
  - corY[]  (2 elements)

# build mcmc
vbgf1d.mcmc <- buildMCMC(vbgf1d.conf)
# compile mcmc and specify the project model
vbgf1d.mcmc.c <- compileNimble(vbgf1d.mcmc)#, project = vbgf1b.c)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Samples

removing output: “warning: logProb of data node length[7]: logProb less than -1e12.”

# sample model
vbgf1d.samp <- runMCMC(vbgf1d.mcmc.c, nchains = 2, niter = 350000, nburnin = 340000, thin=2, samplesAsCodaMCMC = TRUE, WAIC=TRUE)

# vbgf1d.samp <- runMCMC(vbgf1d.mcmc.c, nchains = 2, niter = 25000, nburnin = 20000, thin=2, samplesAsCodaMCMC = TRUE, WAIC=TRUE)

summary(vbgf1d.samp$samples)

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 2 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean        SD  Naive SE Time-series SE
b1[1]         1.23177 9.167e+00 9.167e-02      5.517e-01
b1[2]        -1.79722 9.820e+00 9.820e-02      2.439e-01
b1[3]         4.14670 9.926e+00 9.926e-02      5.987e-01
b1[4]         0.52597 9.948e+00 9.948e-02      1.956e-01
b1[5]        -0.44423 1.009e+01 1.009e-01      1.557e-01
b1[6]         0.32384 1.002e+01 1.002e-01      2.060e-01
b1[7]         0.31509 9.974e+00 9.974e-02      2.728e-01
b1[8]         7.61613 9.904e+00 9.904e-02      9.188e-01
b1[9]        -1.28232 9.579e+00 9.579e-02      5.597e-01
b1[10]        6.17598 1.005e+01 1.005e-01      4.840e-01
b1[11]       -0.06338 9.701e+00 9.701e-02      3.127e-01
b2[1]        -0.04140 3.035e-02 3.035e-04      5.397e-03
b2[2]         0.61062 8.496e-01 8.496e-03      3.241e-01
b2[3]         0.10257 4.994e-03 4.994e-05      5.644e-04
b2[4]         0.01933 3.525e-03 3.525e-05      4.821e-04
b2[5]         0.50818 6.432e-01 6.432e-03      1.139e-01
b2[6]         0.01401 7.895e-02 7.895e-04      9.929e-03
b2[7]         0.02723 2.889e-03 2.889e-05      4.485e-04
b2[8]        -0.15067 3.099e-02 3.099e-04      1.133e-02
b2[9]        -0.17601 3.204e-02 3.204e-04      6.100e-03
b2[10]       -0.01995 2.202e-02 2.202e-04      3.150e-03
b2[11]       -0.06064 4.839e-02 4.839e-04      1.015e-02
b3[1]        -0.32302 1.604e-01 1.604e-03      2.181e-02
b3[2]         1.22675 5.063e-01 5.063e-03      7.370e-02
b3[3]         3.95586 1.966e-01 1.966e-03      5.493e-02
b3[4]         3.45531 4.138e-01 4.138e-03      8.538e-02
b3[5]         1.82345 1.560e+00 1.560e-02      3.108e-01
b3[6]         0.09555 8.023e-01 8.023e-03      1.203e-01
b3[7]         4.02719 3.112e-01 3.112e-03      8.222e-02
b3[8]        -1.69615 2.620e-01 2.620e-03      5.220e-02
b3[9]        -1.37836 3.335e-01 3.335e-03      6.563e-02
b3[10]       -0.26280 2.913e-01 2.913e-03      3.761e-02
b3[11]       -0.70708 4.739e-01 4.739e-03      1.373e-01
klt[1, 1]   955.48481 2.947e+01 2.947e-01      3.374e+00
klt[2, 1]  1039.06675 4.514e+01 4.514e-01      6.132e+00
klt[3, 1]  1630.68482 4.037e+01 4.037e-01      1.032e+01
klt[4, 1]  2379.37317 8.368e+01 8.368e-01      2.674e+01
klt[5, 1]  1053.69578 1.481e+02 1.481e+00      3.408e+01
klt[6, 1]  1166.35657 6.498e+01 6.498e-01      4.450e+00
klt[7, 1]  2303.78578 6.478e+01 6.478e-01      1.122e+01
klt[8, 1]   876.94948 2.488e+01 2.488e-01      8.303e+00
klt[9, 1]   907.65532 1.854e+01 1.854e-01      2.849e+00
klt[10, 1]  982.78143 3.807e+01 3.807e-01      3.230e+00
klt[11, 1] 1034.85224 9.745e+01 9.745e-01      2.553e+01
klt[1, 2]     0.47455 6.526e-02 6.526e-04      7.688e-03
klt[2, 2]     0.65016 1.816e-01 1.816e-03      2.747e-02
klt[3, 2]     0.10338 7.016e-03 7.016e-05      2.262e-03
klt[4, 2]     0.04878 3.127e-03 3.127e-05      9.685e-04
klt[5, 2]     0.43744 2.399e-01 2.399e-03      5.093e-02
klt[6, 2]     0.37816 1.074e-01 1.074e-03      1.034e-02
klt[7, 2]     0.04700 2.329e-03 2.329e-05      4.954e-04
klt[8, 2]     0.41537 5.029e-02 5.029e-04      2.225e-02
klt[9, 2]     0.52400 4.623e-02 4.623e-04      8.997e-03
klt[10, 2]    0.28429 4.313e-02 4.313e-04      4.379e-03
klt[11, 2]    0.43919 1.234e-01 1.234e-03      3.608e-02
klt[1, 3]    -0.49320 2.247e-01 2.247e-03      2.452e-02
klt[2, 3]     0.35797 5.540e-01 5.540e-03      7.435e-02
klt[3, 3]    -2.24078 2.276e-01 2.276e-03      6.826e-02
klt[4, 3]    -4.05277 4.130e-01 4.130e-03      1.067e-01
klt[5, 3]    -0.51270 1.488e+00 1.488e-02      2.794e-01
klt[6, 3]     0.11521 8.249e-01 8.249e-03      1.050e-01
klt[7, 3]    -4.37827 2.770e-01 2.770e-03      1.109e-01
klt[8, 3]    -1.01559 2.525e-01 2.525e-03      8.033e-02
klt[9, 3]    -0.26150 1.291e-01 1.291e-03      1.849e-02
klt[10, 3]   -0.84863 3.552e-01 3.552e-03      3.214e-02
klt[11, 3]   -0.26311 3.207e-01 3.207e-03      7.445e-02

2. Quantiles for each variable:

                 2.5%        25%        50%        75%      97.5%
b1[1]       -15.84518   -5.30811    1.08608  7.642e+00  1.944e+01
b1[2]       -20.49907   -8.52108   -1.73102  4.776e+00  1.794e+01
b1[3]       -14.25120   -2.71836    3.73561  1.084e+01  2.409e+01
b1[4]       -19.31703   -6.18144    0.67607  7.325e+00  1.965e+01
b1[5]       -20.46685   -7.20781   -0.56139  6.398e+00  1.908e+01
b1[6]       -19.13388   -6.52037    0.31828  6.955e+00  2.017e+01
b1[7]       -19.17323   -6.29891    0.18486  6.966e+00  2.005e+01
b1[8]       -10.20372    0.57257    7.32230  1.425e+01  2.783e+01
b1[9]       -20.09179   -7.92539   -1.11219  5.237e+00  1.723e+01
b1[10]      -13.44814   -0.69656    6.05705  1.293e+01  2.587e+01
b1[11]      -18.43861   -6.70025   -0.18296  6.475e+00  1.902e+01
b2[1]        -0.10081   -0.06355   -0.03929 -2.047e-02  1.916e-02
b2[2]         0.04238    0.25718    0.39385  5.922e-01  4.007e+00
b2[3]         0.09452    0.09918    0.10191  1.051e-01  1.152e-01
b2[4]         0.01256    0.01701    0.01929  2.138e-02  2.701e-02
b2[5]        -0.15176    0.11055    0.28733  6.652e-01  2.497e+00
b2[6]        -0.15004   -0.03585    0.01843  6.917e-02  1.561e-01
b2[7]         0.02228    0.02503    0.02704  2.888e-02  3.340e-02
b2[8]        -0.20302   -0.17100   -0.15558 -1.304e-01 -9.105e-02
b2[9]        -0.23961   -0.19761   -0.17467 -1.530e-01 -1.147e-01
b2[10]       -0.06789   -0.03215   -0.01810 -4.933e-03  1.588e-02
b2[11]       -0.15252   -0.09654   -0.06032 -2.740e-02  3.836e-02
b3[1]        -0.62287   -0.43665   -0.32780 -2.177e-01  4.966e-03
b3[2]         0.24224    0.88972    1.22994  1.596e+00  2.160e+00
b3[3]         3.51275    3.81678    4.01145  4.088e+00  4.249e+00
b3[4]         2.60956    3.19782    3.48566  3.755e+00  4.208e+00
b3[5]        -1.24013    0.81526    1.62145  2.787e+00  5.582e+00
b3[6]        -1.32410   -0.48514    0.03103  6.719e-01  1.716e+00
b3[7]         3.45136    3.79124    4.01879  4.275e+00  4.582e+00
b3[8]        -2.22000   -1.85165   -1.69868 -1.542e+00 -1.164e+00
b3[9]        -2.19732   -1.55779   -1.34616 -1.150e+00 -8.218e-01
b3[10]       -0.88196   -0.45318   -0.25031 -6.288e-02  2.571e-01
b3[11]       -1.93837   -0.95182   -0.61409 -3.544e-01 -3.119e-02
klt[1, 1]   909.67046  933.96471  951.37115  9.721e+02  1.025e+03
klt[2, 1]   970.30300 1007.60103 1032.51643  1.061e+03  1.150e+03
klt[3, 1]  1534.37232 1597.94737 1655.95381  1.661e+03  1.665e+03
klt[4, 1]  2247.38352 2319.53963 2371.64916  2.429e+03  2.561e+03
klt[5, 1]   886.59641  950.45294 1009.91866  1.109e+03  1.487e+03
klt[6, 1]  1063.28607 1119.36460 1158.33067  1.204e+03  1.318e+03
klt[7, 1]  2167.50113 2248.31078 2358.72050  2.359e+03  2.359e+03
klt[8, 1]   845.90233  859.63656  871.80564  8.840e+02  9.409e+02
klt[9, 1]   874.02153  895.06497  906.98867  9.197e+02  9.463e+02
klt[10, 1]  918.34284  955.81385  977.76120  1.006e+03  1.066e+03
klt[11, 1]  904.52432  950.11929 1015.24340  1.109e+03  1.225e+03
klt[1, 2]     0.34642    0.42897    0.47443  5.208e-01  5.979e-01
klt[2, 2]     0.32653    0.52464    0.63167  7.708e-01  1.021e+00
klt[3, 2]     0.09606    0.09825    0.09957  1.080e-01  1.188e-01
klt[4, 2]     0.04287    0.04643    0.04898  5.128e-02  5.415e-02
klt[5, 2]     0.12251    0.25198    0.37072  5.806e-01  1.022e+00
klt[6, 2]     0.21174    0.29140    0.36722  4.492e-01  6.163e-01
klt[7, 2]     0.04368    0.04512    0.04644  4.891e-02  5.141e-02
klt[8, 2]     0.29943    0.39461    0.41924  4.517e-01  4.888e-01
klt[9, 2]     0.43832    0.49187    0.52374  5.515e-01  6.172e-01
klt[10, 2]    0.20766    0.25222    0.28337  3.124e-01  3.783e-01
klt[11, 2]    0.25909    0.32762    0.42755  5.438e-01  6.581e-01
klt[1, 3]    -0.97031   -0.64327   -0.47252 -3.232e-01 -1.136e-01
klt[2, 3]    -0.98525    0.03237    0.43367  7.853e-01  1.182e+00
klt[3, 3]    -2.56362   -2.40161   -2.31561 -2.091e+00 -1.761e+00
klt[4, 3]    -4.71562   -4.39626   -4.10546 -3.738e+00 -3.240e+00
klt[5, 3]    -3.68564   -1.46893   -0.44674  5.924e-01  1.889e+00
klt[6, 3]    -1.65420   -0.44982    0.22469  7.547e-01  1.433e+00
klt[7, 3]    -4.80539   -4.65113   -4.36529 -4.187e+00 -3.819e+00
klt[8, 3]    -1.66254   -1.10496   -0.97216 -8.353e-01 -6.886e-01
klt[9, 3]    -0.54240   -0.34199   -0.25214 -1.744e-01 -2.686e-02
klt[10, 3]   -1.51632   -1.10563   -0.83432 -5.954e-01 -1.724e-01
klt[11, 3]   -0.85166   -0.53630   -0.24406  1.653e-02  2.574e-01

mcmc_trace(vbgf1d.samp$samples)

gelman.diag(vbgf1d.samp$samples)

Potential scale reduction factors:

           Point est. Upper C.I.
b1[1]            1.02       1.05
b1[2]            1.00       1.02
b1[3]            1.04       1.18
b1[4]            1.00       1.01
b1[5]            1.00       1.00
b1[6]            1.00       1.00
b1[7]            1.01       1.03
b1[8]            1.01       1.03
b1[9]            1.08       1.32
b1[10]           1.00       1.00
b1[11]           1.04       1.15
b2[1]            1.11       1.27
b2[2]            1.05       1.19
b2[3]            1.62       2.89
b2[4]            1.02       1.07
b2[5]            1.36       2.25
b2[6]            1.09       1.11
b2[7]            1.03       1.14
b2[8]            1.10       1.33
b2[9]            1.41       2.38
b2[10]           1.00       1.01
b2[11]           1.17       1.57
b3[1]            1.04       1.09
b3[2]            1.05       1.19
b3[3]            2.43       6.89
b3[4]            1.07       1.21
b3[5]            1.25       1.81
b3[6]            1.06       1.07
b3[7]            2.95       6.68
b3[8]            1.06       1.21
b3[9]            1.53       3.06
b3[10]           1.01       1.02
b3[11]           1.46       2.37
klt[1, 1]        1.04       1.08
klt[2, 1]        1.07       1.17
klt[3, 1]        4.69      26.19
klt[4, 1]        1.28       2.27
klt[5, 1]        1.28       2.23
klt[6, 1]        1.05       1.19
klt[7, 1]        3.91      22.79
klt[8, 1]        1.04       1.18
klt[9, 1]        1.09       1.27
klt[10, 1]       1.07       1.30
klt[11, 1]       1.52       3.12
klt[1, 2]        1.05       1.15
klt[2, 2]        1.05       1.06
klt[3, 2]        3.79      17.27
klt[4, 2]        1.59       2.93
klt[5, 2]        1.08       1.10
klt[6, 2]        1.04       1.12
klt[7, 2]        4.16      18.46
klt[8, 2]        1.06       1.25
klt[9, 2]        1.05       1.10
klt[10, 2]       1.05       1.20
klt[11, 2]       1.30       2.01
klt[1, 3]        1.04       1.11
klt[2, 3]        1.02       1.04
klt[3, 3]        2.81       8.48
klt[4, 3]        1.22       1.81
klt[5, 3]        1.03       1.12
klt[6, 3]        1.02       1.02
klt[7, 3]        3.72      11.61
klt[8, 3]        1.08       1.31
klt[9, 3]        1.02       1.03
klt[10, 3]       1.02       1.09
klt[11, 3]       1.33       2.16

Multivariate psrf

10.2

effectiveSize(vbgf1d.samp$samples)

      b1[1]       b1[2]       b1[3]       b1[4]       b1[5]       b1[6] 
 277.303570 1655.795048  283.940779 2593.234206 4202.864253 2372.884782 
      b1[7]       b1[8]       b1[9]      b1[10]      b1[11]       b2[1] 
1342.982141  114.042927  279.825441  483.400284  953.575140   31.778848 
      b2[2]       b2[3]       b2[4]       b2[5]       b2[6]       b2[7] 
  58.742815   82.528124   52.688560   46.075704   69.487666   38.421280 
      b2[8]       b2[9]      b2[10]      b2[11]       b3[1]       b3[2] 
  10.855207   23.331483   50.736635   22.425632   54.300314   45.506735 
      b3[3]       b3[4]       b3[5]       b3[6]       b3[7]       b3[8] 
  18.871407   22.527719   23.290540   49.682297   14.813583   26.472603 
      b3[9]      b3[10]      b3[11]   klt[1, 1]   klt[2, 1]   klt[3, 1] 
  24.616606   60.736811   16.148245   76.347788   50.514005    6.064103 
  klt[4, 1]   klt[5, 1]   klt[6, 1]   klt[7, 1]   klt[8, 1]   klt[9, 1] 
   8.280631   45.419129  214.896662    5.381439   17.695268   51.547888 
 klt[10, 1]  klt[11, 1]   klt[1, 2]   klt[2, 2]   klt[3, 2]   klt[4, 2] 
 144.435366   14.892376   76.334946   56.131392   18.087843   22.390441 
  klt[5, 2]   klt[6, 2]   klt[7, 2]   klt[8, 2]   klt[9, 2]  klt[10, 2] 
  39.834436  109.944925    8.736615   10.560269   35.657337   99.288503 
 klt[11, 2]   klt[1, 3]   klt[2, 3]   klt[3, 3]   klt[4, 3]   klt[5, 3] 
  14.432407   85.141868   52.576788   17.443881   23.840913   24.875412 
  klt[6, 3]   klt[7, 3]   klt[8, 3]   klt[9, 3]  klt[10, 3]  klt[11, 3] 
  66.854249    7.155436   21.350042   60.571555  123.632214   20.919353 

# vbgf1d.samp$samples %>%
#   gather_draws(klt[spat.unit,par], sep = ",") %>%
#   filter(par == 1) %>%
#   ggplot() + 
#   geom_density(aes(x = .value, color = .variable)) + # age parameter i
#   facet_wrap(~spat.unit, scales = "free", nrow = 3) +
#   theme_light() +
#   
# vbgf1d.samp$samples %>%
#   gather_draws(klt[spat.unit,par], b1[spat.unit], sep = ",") %>%
#   filter(par == 1 | is.na(par)) %>%
#   ggplot() + 
#   geom_density(aes(x = .value, color = .variable)) + # age parameter i
#   facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
#   theme_light() 
#   
# vbgf1d.samp$samples %>%
#   gather_draws(klt[spat.unit,par], b1[spat.unit], b2[spat.unit], b3[spat.unit], sep = ",") %>%
#   filter(par == 2) %>%
#   ggplot() + 
#   geom_density(aes(x = .value, color = factor(par))) + # age parameter i
#   facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
#   theme_light() +
#   
# vbgf1d.samp$samples %>%
#   gather_draws(klt[spat.unit,par], beta1, sep = ",") %>%
#   filter(par == 3) %>%
#   ggplot() + 
#   geom_density(aes(x = .value, color = factor(par))) + # age parameter i
#   facet_wrap(par~spat.unit, scales = "free", nrow = 3) +
#   theme_light() 

e. Compare WAIC with mnorm and with/without t0.

c(c = vbgf1c.samp$WAIC, d = vbgf1d.samp$WAIC)

$c
nimbleList object of type waicNimbleList
Field "WAIC":
[1] 309008.4
Field "lppd":
[1] -154463.6
Field "pWAIC":
[1] 40.60854

$d
nimbleList object of type waicNimbleList
Field "WAIC":
[1] 308291.9
Field "lppd":
[1] -154094.3
Field "pWAIC":
[1] 51.64428